facs.summary.plot<-function(x, ncells=10000, view.plot=T){

  ##tick_locations = c(0, 128, 192, 256, 384, 448, 512, 640, 704, 768, 896, 960, 1024)  
  ##tick_locations = c(0, 128, 256, 384, 512, 640, 768, 896, 1024)  
  ##tick_locations = c(0, 187.9363, 256, 434.9363, 512, 690.9363, 768, 946.9363, 1024)  
  tick_locations = c(0, 187.9363, 256, 434.9363, 512, 768, 1024)  
  tick_labels = round(10^(tick_locations/256))
  
  z<-list() ## for saving results
  xx<-exprs(x)  
  ## delete false
  del.idx<-which(xx[,"SSC-H"]<=0 | xx[,"FL1-H"]<=0 | xx[,"FSC-H"]<=0)
  ##del.idx<-c(del.idx, which())
  if(length(del.idx)>0){
    xx<-xx[-del.idx,]
  }
  
  
  ## select 10000 cells from
  ##ncells<-15000
  ## select the most dense area
  fsc.density<-density(xx[,"FSC-H"])
  fsc.max<-fsc.density$x[which.max(fsc.density$y)]
  ##fsc.max<-256
  ssc.density<-density(xx[,"SSC-H"])
  ssc.max<-ssc.density$x[which.max(ssc.density$y)]
  ##fl1 density
  fl1.density<-density(xx[,"FL1-H"])
  fl1.max<-fl1.density$x[which.max(fl1.density$y)]
  ##ssc.max<-356
  
  d<-sqrt((xx[,"FSC-H"]-fsc.max)^2+(xx[,"SSC-H"]-ssc.max)^2)
  o<-order(d)
  if(length(o)>ncells){
    sel.idx<-o[1:ncells]
  }else{
    sel.idx<-o[1:length(o)]
  }
  
  sel.xx<-xx[sel.idx,]
  ## fsc
  sel.fsc.density<-density(sel.xx[,"FSC-H"])
  sel.fsc.max<-sel.fsc.density$x[which.max(sel.fsc.density$y)]
  ## ssc
  sel.ssc.density<-density(sel.xx[,"SSC-H"])
  sel.ssc.max<-sel.ssc.density$x[which.max(sel.ssc.density$y)]
  ## fl1
  sel.fl1.density<-density(sel.xx[,"FL1-H"])
  sel.fl1.max<-sel.fl1.density$x[which.max(sel.fl1.density$y)]
  
  
  ##cat(fsc.max, ssc.max, "mean:", mean(sel.xx[,"FL1-H"]), "dim:", dim(sel.xx), "\n");flush.console()
  cat("Total.cells:", nrow(xx), "\n");flush.console()
  
  z$fsc.max<-fsc.max
  z$ssc.max<-ssc.max
  z$sel.fsc.max<-sel.fsc.max
  z$sel.ssc.max<-sel.ssc.max
  z$total.cells<-nrow(xx)
  
  ### ----- plot ------
  if(view.plot==T){
    ##par(mfrow=c(3,2))
    
    layout(matrix(c(1:10), 5, 2, byrow = FALSE))
    
    ##layout(matrix(c(1:2), 2, 1, byrow = FALSE))
    
    ## FSC vs. SSC and the automatically selected 10000 cells
    plot(xx[,"FSC-H"], xx[,"SSC-H"], pch=20, col="#55555511", xlab="FSC", ylab="SSC", main="Scatter plot (FSC vs. SSC)", xaxt = "n", yaxt="n", xlim=c(0,1024), ylim=c(0,1024), bty="n")
    axis(1, at = tick_locations, labels = tick_labels)
    axis(2, at = tick_locations, labels = tick_labels)
    ##points(sel.xx[,"FSC-H"], sel.xx[,"SSC-H"], pch=19, col="#0000ff03", xlab="FSC", ylab="SSC", main="Scatter plot (FSC vs. SSC)", xlim=c(0,500), ylim=c(0,800))
    
    ## FSC vs. FL1 without the selected 10000 cells
    plot(xx[,"FSC-H"], xx[,"FL1-H"], pch=20, col="#00005511", xlab="FSC", ylab="FL1", main="Scatter plot (FSC vs. FL1)", xaxt = "n", yaxt="n", xlim=c(0,1024), ylim=c(0,1024), bty="n")
    axis(1, at = tick_locations, labels = tick_labels)
    axis(2, at = tick_locations, labels = tick_labels)
    ##points(xx[,"FSC-H"], xx[,"FL1-H"], pch=19, col="#0000ff", xlab="FSC", ylab="FL1", main="Scatter plot (FSC vs. SSC)", xlim=c(0,500), ylim=c(0,800))
    
    
    ## FSC vs cell density
    plot(density(xx[,"FSC-H"]), col="#000000", ylab="density", xlab="FSC", xaxt="n", main="FSC density", xlim=c(0,1024), bty="n")
    axis(1, at = tick_locations, labels = tick_labels)
    ##plot(density(log(10^(xx[,"FL1-H"]/256)), br=100))
    abline(v=fsc.max, lty=2)
    text(fsc.max, max(fsc.density$y)*0.1, round(10^(fsc.max/256),2), cex=2, pos=4)
    
    ## SSC vs. cell density
    plot(density(xx[,"SSC-H"]), col="#000000", ylab="density", xlab="SSC", xaxt="n", main="SSC density", xlim=c(0,1024), bty="n")
    axis(1, at = tick_locations, labels = tick_labels)
    abline(v=ssc.max, lty=2)
    text(ssc.max, max(ssc.density$y)*0.1, round(10^(ssc.max/256),2), cex=2, pos=4)

    ### FL1 of the selected cells vs. their density
    plot(density(xx[,"FL1-H"]), col="#000000", ylab="density", xlab="FL1", bty="n", main="FL1 density", xaxt="n", xlim=c(0,1024), bty="n")
    axis(1, at = tick_locations, labels = tick_labels)
    abline(v=fl1.max, lty=2)
    text(fl1.max, max(fl1.density$y)*0.1, round(10^(fl1.max/256),2), cex=2, pos=4)
    

    ##======= 1st selected
    
    
    ## FSC vs. SSC and the automatically selected 10000 cells
    plot(xx[,"FSC-H"], xx[,"SSC-H"], pch=20, col="#55555511", xlab="FSC", ylab="SSC", main="Scatter plot (FSC vs. SSC)", xaxt = "n", yaxt="n", xlim=c(0,1024), ylim=c(0,1024), bty="n")
    axis(1, at = tick_locations, labels = tick_labels)
    axis(2, at = tick_locations, labels = tick_labels)
    points(sel.xx[,"FSC-H"], sel.xx[,"SSC-H"], pch=19, col="#ff000011", xlab="FSC", ylab="SSC", main="Scatter plot (FSC vs. SSC)", xlim=c(0,500), ylim=c(0,800))
    
    ## FSC vs. FL1 without the selected 10000 cells
    plot(xx[,"FSC-H"], xx[,"FL1-H"], pch=20, col="#55555511", xlab="FSC", ylab="FL1", main="Scatter plot (FSC vs. FL1)", xaxt = "n", yaxt="n", xlim=c(0,1024), ylim=c(0,1024), bty="n")
    axis(1, at = tick_locations, labels = tick_labels)
    axis(2, at = tick_locations, labels = tick_labels)
    points(sel.xx[,"FSC-H"], sel.xx[,"FL1-H"], pch=19, col="#ff000011", xlab="FSC", ylab="FL1", main="Scatter plot (FSC vs. SSC)", xlim=c(0,500), ylim=c(0,800))
    
    
    ## FSC vs cell density
    plot(density(sel.xx[,"FSC-H"]), col="#000000", ylab="density", xlab="FSC", xaxt="n", main="FSC density", xlim=c(0,1024), bty="n")
    axis(1, at = tick_locations, labels = tick_labels)
    ##plot(density(log(10^(xx[,"FL1-H"]/256)), br=100))
    abline(v=sel.fsc.max, lty=2)
    text(sel.fsc.max, max(sel.fsc.density$y)*0.1, round(10^(sel.fsc.max/256),2), cex=2, pos=4)
    
    ## SSC vs. cell density
    plot(density(sel.xx[,"SSC-H"]), col="#000000", ylab="density", xlab="SSC", xaxt="n", main="SSC density", xlim=c(0,1024), bty="n")
    axis(1, at = tick_locations, labels = tick_labels)
    abline(v=sel.ssc.max, lty=2)
    text(sel.ssc.max, max(sel.ssc.density$y)*0.1, round(10^(sel.ssc.max/256),2), cex=2, pos=4)
    
    ### FL1 of the selected cells vs. their density
    plot(density(sel.xx[,"FL1-H"]), col="#000000", ylab="density", xlab="FL1", bty="n", main="FL1 density", xaxt="n", xlim=c(0,1024), bty="n")
    axis(1, at = tick_locations, labels = tick_labels)
    abline(v=sel.fl1.max, lty=2)
    text(sel.fl1.max, max(sel.fl1.density$y)*0.1, round(10^(sel.fl1.max/256),2), cex=2, pos=4)
    
  }
  
  z$sel.xx<-sel.xx
  return(z)
  
}




facs.compare.plot<-function(zz, ncells=10000, plot.type="FSC", overlap=FALSE, selection=TRUE){
  
  tick_locations = c(0, 187.9363, 256, 434.9363, 512, 768, 1024)  
  tick_labels = round(10^(tick_locations/256))
  
  nplot<-length(zz)
  sample.names<-names(zz)
  z<-list() ## for saving results
  x11(width=5, height=17)
  layout(matrix(c(1:nplot), nplot, 1, byrow = FALSE))
  
  for(i in 1:nplot){
    x<-zz[[i]]
    xx<-exprs(x)  
    ## delete false 
    del.idx<-which(xx[,"SSC-H"]<=0 | xx[,"FL1-H"]<=0 | xx[,"FSC-H"]<=0)
    ##del.idx<-c(del.idx, which())
    if(length(del.idx)>0){
      xx<-xx[-del.idx,]
    }
    
    if(selection==TRUE){
      ##ncells<-15000
      ## select the most dense area
      fsc.density<-density(xx[,"FSC-H"])
      fsc.max<-fsc.density$x[which.max(fsc.density$y)]
      ##fsc.max<-256
      ssc.density<-density(xx[,"SSC-H"])
      ssc.max<-ssc.density$x[which.max(ssc.density$y)]
      ##fl1 density
      fl1.density<-density(xx[,"FL1-H"])
      fl1.max<-fl1.density$x[which.max(fl1.density$y)]
      d<-sqrt((xx[,"FSC-H"]-fsc.max)^2+(xx[,"SSC-H"]-ssc.max)^2)
      o<-order(d)
      if(length(o)>ncells){
        sel.idx<-o[1:ncells]
      }else{
        sel.idx<-o[1:length(o)]
      }
      sel.xx<-xx[sel.idx,]
    }else{
      sel.xx<-xx
    }
    
    tt<-paste(plot.type, "-H", sep="")
    sel.density<-density(sel.xx[,tt])
    sel.val.max<-sel.density$x[which.max(sel.density$y)]
    
    cat("Number of cells:", nrow(sel.xx), "\n");flush.console()
    ### ----- plot ------
    
    plot(density(sel.xx[,tt]), col="#000000", ylab="density", xlab=plot.type, xaxt="n", main=paste(tt, sample.names[i]), xlim=c(0,1024), bty="n")
    axis(1, at = tick_locations, labels = tick_labels)
    abline(v=sel.val.max, lty=2)
    text(sel.val.max, max(sel.density$y)*0.1, round(10^(sel.val.max/256),2), cex=2, pos=4)
  }
    
  ##z$sel.xx<-sel.xx
  return(z)
  
}



facs.sumcom.plot<-function(zz, ncells=10000, plot.type="FSC", overlap=FALSE, selection=TRUE){
  
  tick_locations = as.numeric(names(zz))
  tick_labels = as.numeric(names(zz)) 
  ##round(10^(tick_locations/256))
  
  nplot<-length(zz)
  sample.names<-names(zz)
  z<-list() ## for saving results
  x11(width=7, height=14)
  layout(matrix(c(1:2), 2, 1, byrow = FALSE))
  
  max.arr<-c()
  var.arr<-c()
  for(i in 1:nplot){
    x<-zz[[i]]
    xx<-exprs(x)  
    ## delete false 
    del.idx<-which(xx[,"SSC-H"]<=0 | xx[,"FL1-H"]<=0 | xx[,"FSC-H"]<=0)
    ##del.idx<-c(del.idx, which())
    if(length(del.idx)>0){
      xx<-xx[-del.idx,]
    }
    
    if(selection==TRUE){
      ##ncells<-15000
      ## select the most dense area
      fsc.density<-density(xx[,"FSC-H"])
      fsc.max<-fsc.density$x[which.max(fsc.density$y)]
      ##fsc.max<-256
      ssc.density<-density(xx[,"SSC-H"])
      ssc.max<-ssc.density$x[which.max(ssc.density$y)]
      ##fl1 density
      fl1.density<-density(xx[,"FL1-H"])
      fl1.max<-fl1.density$x[which.max(fl1.density$y)]
      d<-sqrt((xx[,"FSC-H"]-fsc.max)^2+(xx[,"SSC-H"]-ssc.max)^2)
      o<-order(d)
      if(length(o)>ncells){
        sel.idx<-o[1:ncells]
      }else{
        sel.idx<-o[1:length(o)]
      }
      sel.xx<-xx[sel.idx,]
    }else{
      sel.xx<-xx
    }
    
    tt<-paste(plot.type, "-H", sep="")
    sel.density<-density(sel.xx[,tt])
    sel.val.max<-sel.density$x[which.max(sel.density$y)]
    
    cat("Number of cells:", nrow(sel.xx), "\n");flush.console()
    ### ----- plot ------
    max.arr<-c(max.arr, sel.val.max)
    var.arr<-c(var.arr, var(sel.xx[,tt]))
    ##plot(density(sel.xx[,tt]), col="#000000", ylab="density", xlab=plot.type, xaxt="n", main=paste(tt, sample.names[i]), xlim=c(0,1024), bty="n")
    ##axis(1, at = tick_locations, labels = tick_labels)
    ##abline(v=sel.val.max, lty=2)
    ##text(sel.val.max, max(sel.density$y)*0.1, round(10^(sel.val.max/256),2), cex=2, pos=4)

  }
  z$max.arr<-max.arr
  z$var.arr<-var.arr
  plot(x=tick_locations, y=max.arr, ylim=c(min(max.arr)*0.8, max(max.arr)*1.2), xaxt="n", bty="n")
  axis(1, at = tick_locations, labels = tick_labels)
  plot(x=tick_locations, y=var.arr, ylim=c(min(var.arr)*0.8, max(var.arr)*1.2), xaxt="n", bty="n")
  axis(1, at = tick_locations, labels = tick_labels)
  return(z)
  
}

